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SUMMARY 

Very  large  sets  of  simultaneous  linear  equations  can  usually  be 
solved  to  reasonable  accuracy  more  quickly  by  an  iterative  method  than  by 
a  direct  one,  but  most  iterative  methods  only  converge  if  the  elements  of 
the  principal  diagonal  of  the  coefficient  matrix  dominate  all  the  other 
elements.  This  paper  describes  an  iterative  method  which  depends  for  its 
convergence  only  on  dominance  by  the  elements  of  an  arbitrarily  wide  band 
of  diagonals  centred  on  the  principal  diagonal.  The  theory  of  the  method 
is  fully  worked  out,  and  a  computer  program  implementing  it  in  Fortran  IV 
is  presented. 
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1  INTRODUCTION 

Boundary  value  problems  in  field  theory,  including  those  of  the 
electromagnetic  radiation  from  antennas,  are  conveniently  expressible  in 
4  the  form  of  linear  integral  equations.  It  has  become  customary  to  solve 

these  equations  using  the  so-called  "method  of  moments"  (Ref  1) .  This 

*  method  consists  of  expanding  the  unknown  field  function  as  a  linear 
combination  of  suitably  chosen  functions  and  taking  the  integral  average  of 
the  resulting  form  of  the  integral  equation  with  each  member  of  the  set  of 
functions  (or  sometimes  with  each  member  of  another  set  of  functions).  The 
integral  equation  then  reduces  to  a  fully-determined  set  of  linear 
simultaneous  equations  in  as  many  unknown  coefficients  as  there  are  functions 
in  the  set  (if  2  sets  of  functions  are  used  there  should  be  the  same  number 
of  functions  in  each  set) . 

In  principle,  the  solution  of  a  fully-determined  set  of  linear 

simultaneous  equations  is  an  elementary  problem.  It  is  known,  however,  that 

the  number  of  arithmetic  operations  required  to  solve  exactly  N  linear 

3 

simultaneous  equations  in  N  unknowns  is  proportional  to  N  (for  large  N), 
and  if  N  is  allowed  to  increase  substantially  beyond  about  100  the  solution 

*  time  rapidly  becomes  intolerably  long  except  on  very  fast  computers. 

Moreover,  the  accuracy  of  such  a  calculation  is  limited  only  by  the 

%  precision  with  which  the  coefficients  are  represented,  so  that  the  solution 


1 


may  be  given  to  10  or  20  significant  figures  although  all  but  the  first  3 
or  4  of  these  are  invalidated  by  the  inaccuracy  of  the  physical  model  and 
the  mathematical  approximations  employed. 

The  use  of  an  iterative  method  of  solution  addresses  both  these 

objections.  As  will  be  shown  later,  the  number  of  arithmetic  operations 

required  to  perform  one  iteration  on  an  N  x  N  matrix  of  coefficients  is 
2 

approximately  N  ,  and  so  the  solution  time  can  in  principle  be  substantially 
reduced  if  the  number  of  iterations  required  is  at  least  several  times 
smaller  than  the  order  N  of  the  matrix.  This  is  generally  possible  whenever 
the  solution  is  not  required  to  attain  great  accuracy,  which  is  precisely 
the  case  here.  This  paper  describes  an  iterative  method  using  banded 
matrices  which,  unlike  the  methods  commonly  used,  does  not  rely  on  diagonal 
dominance  in  the  coefficient  matrix,  but  only  on  dominance  by  an  arbitrarily 
wide  band  of  diagonals  centred  on  the  principal  diagonal. 


2  THE  ITERATIVE  PRINCIPLE 


The  application  of  iteration  to  the  problem  of  solving  N  linear 
simultaneous  equations  in  N  unknowns  is  most  conveniently  described  in  the 
notation  of  matrix  algebra.  Let  the  system  be  represented  by  the  matrix 
equation 

Ax  =  f  (1) 


where  A  is  an  N  x  N  matrix  of  coefficients  and  x  and  f  are  N  x  1  matrices 

(N  -  dimensional  column  vectors)  representing  the  unknowns  and  the  right- 

hand  sides  of  the  equations,  respectively.  The  iterative  principle  is 

applied  to  this  system  by  separating  the  matrix  A  into  2  parts,  one 

containing  most  of  the  large  elements  of  A  so  that  it  will  dominate  the  other 

part.  If  we  refer  to  the  parts  as  the  large  and  the  small  part,  denoted  by 

A,  and  A  ,  we  then  have 
1  s 

A  =  Aj  +  A^  or  (A^  +  Ag)  x  =  f  from  equation  (1),  so  (2) 

we  may  iterate  using  the  equation 


V(i>  =  f  -  A  x(i_1) 
1  s 


(3) 


where  denotes  the  i*"*1  approximation  to  the  exact  solution  x  of  equation 

(1). 
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We  begin  Che  mathematical  analysis  of  this  process  by  deriving  the 
condition  for  it  to  converge.  Let  the  error  of  the  ith  iteration  be  denoted 
by  the  vector  e^,  thus 


e 


(i) 


Then  we  have 


(4) 


=  x(i>  -  x  =  1  (AjX^  -  AjX)  »  Aj  1  (f  -  Asx^  ^  ~  A^x) 

from  equation  (3) 

=  A.  ^((A,  +  A  )x  -  A  x^  ^  -  A,x) 

X  1  s  s  1 

from  equation  (2) 

=  A  *(A  (x  -  x^1  1^)) 

1  s 


(i-1) 


from  equation  (4) ,  and  so 


■<l)  -  (-  A,'1*. 


(5) 


It  then  follows  from  standard  results  in  matrix  theory  (eg  Ref  2,  pp  111-112) 
that  the  necessary  and  sufficient  condition  for  e'1^  to  tend  to  the  null 
vector  as  i  tends  to  infinity  is  that  the  moduli  of  the  latent  roots  of  the 
matrix  (-  A^  ^Ag)  should  all  be  less  than  unity,  and  that  the  convergence  is 
ultimately  of  the  type  normally  called  linear  (the  number  of  correct 
significant  figures  in  x^  increases  linearly  with  i)  since  and 

e^  ^  satisfy  asymptotically  the  equation 


e 


(i) 


Xe 


(i-1) 


(6) 


where  the  scalar  X  is  the  latent  root  largest  in  modulus  of  the  matrix 

(-  A,  ^A  ) .  This  dependence  on  the  latent  roots  of  the  matrix  (-  A,  *A  ) 
is  is 

shows  why  the  parts  of  A  should  be  "large"  and  "small"  to  bring  about 
convergence.  We  can  also  see  from  equation  (5)  that  the  convergence  or 
otherwise  of  the  sequence  of  vectors  e^  does  not  depend  on  so  if  an 

iterative  process  of  the  kind  represented  by  equation  (2)  converges  at  all 
it  will  do  so  from  any  starting  point. 
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BANDED  MATRICES 


If  the  iterative  method  is  to  be  practical,  equation  (3)  must  be 
economically  soluble  at  each  iteration,  so  the  implied  operation  of 
inverting  A^  must  be  arithmetically  much  quicker  than  the  implied  operation 
of  inverting  A  to  solve  equation  (1)  directly.  This  is  commonly  ensured 
by  making  A  triangular  or  diagonal,  as  in  the  Gauss-Seidel  and  Jacobi  methods 
(see  Reference  2,  p  179  et  seq,  for  a  description  of  these),  and  the 
algorithm  for  solving  equation  (3)  is  then  a  very  simple  one.  Unfortunately , 
both  these  methods  diverge  in  general  unless  the  matrix  A  is  diagonally 
dominant,  and,  although  matrices  arising  out  of  antenna  modelling  problems 
usually  have  their  largest  elements  on  their  principal  diagonals,  the  other 
elements  are  not  usually  sufficiently  relatively  small. 

The  technique  described  in  this  paper  addresses  this  problem  by  making 
the  matrix  A^  a  banded  matrix;  that  is,  a  matrix  which  has  null  elements 
everywhere  except  on  the  principal  diagonal  and  on  a  number  of  neighbouring 
sub-diagonals  and  super-diagonals.  (This  idea  is  presented  in  the  paper 
cited  as  Ref  3,  on  which  this  work  is  based.)  It  is  known  (Ref  2,  pp  17-18) 
that  any  square  matrix  may  be  resolved  into  a  product  of  a  lower-triangular 
matrix  and  an  upper-triangular  one,  provided  that  the  square  matrix  and  all 
its  principal  sub-matrices  are  non-singular.  In  antenna  modelling  it  is 
generally  taken  as  an  article  of  faith  that  the  matrix  of  coefficients 
corresponding  to  an  antenna  will  be  non-singular  if  the  exciting  frequency 
does  not  correspond  to  any  of  the  antenna's  free  resonances,  and  that  at  the 
free  resonances  no  unique  solution  is  possible  anyway  (ignoring  the 
difference  between  the  free  resonances  of  the  antenna  and  those  of  the  model 
of  the  antenna).  It  is  then  argued  that,  since  all  the  principal  sub-matrices 
are  the  matrices  associated  with  "sub-problems"  (problems  obtained  by 
omitting  some  of  the  parts  of  the  model),  they  will  also  be  non-singular;  and 
that  the  banded  matrix  A^  and  its  principal  sub-matrices  are  non-singular, 
since  they  arc  obtained  from  the  full  matrix  A  and  its  principal  sub-matrices 
by  simply  ignoring  some  of  the  smaller  interactions  between  parts  of  the 
model.  This  argument  is  patently  unreliable,  but  it  appears  to  be  justified 
in  practice.  Now,  when  the  banded  matrix  A^  is  so  chosen  that  the  number  of 
sub-diagonals  is  equal  to  the  number  of  super-diagonals  (the  common  value 
being  M,  say),  it  turns  out,  as  remarked  in  Ref  3,  that  the  triangular 
matrices  into  which  Aj  may  be  resolved  are  also  banded  and  have  the  same 
bandwidth;  that  is,  the  lower-triangular  matrix  lias  nulls  everywhere  but  on 
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the  principal  diagonal  and  on  the  H  neighbouring  sub-diagonals,  and  the 
upper-triangular  matrix  has  nulls  everywhere  but  on  the  principal  diagonal 
and  on  the  M  neighbouring  super-diagonals.  As  the  matrix  of  coefficients  A 
is  roughly  synmetrical  (in  a  z-j °d  model,  the  action  of  part  i  on  part  j  will 
be  approximately  equal  to  that  of  part  j  on  part  i,  on  account  of  reciprocity) 
it  is  acceptable  and,  indeed,  logical  to  choose  A^  such  that  the  numbers  of 
sub-diagonals  and  super-diagonals  are  equal.  Consequently,  resolving  A^ 
into  a  product  of  banded  triangular  matrices  enables  equation  (3)  to  be 
solved  quickly  at  each  iteration;  the  resolution  only  needs  to  be  done  once 
and  is  arithmetically  much  quicker  than  resolution  of  a  full  matrix  of  the 
same  order  (as  will  be  shown) ,  and  the  solution  of  equation  (3)  becomes  then 
merely  a  matter  of  forward  and  backward  substitution,  made  even  quicker  than 
usual  by  the  presence  of  extra  nulls  in  the  triangular  matrices. 


4  THE  TRIANGULAR  DECOMPOSITION 

The  discussion  in  reference  2  (pp  17-22)  on  triangular  decomposition 

shows  that  it  has  a  certain  arbitrariness;  there  is,  essentially,  one  degree 

of  freedom  in  each  row  (or  each  column)  of  one  of  the  triangular  matrices , 

or  the  degrees  of  freedom  may  be  shared  between  them.  This  is  a  consequence 

2 

of  the  fact  that  the  2  triangular  matrices  between  them  contain  (N  +  N) 

2 

non-null  elements  constrained  only  by  the  (N  )  elements  in  the  original  N  x  N 
matrix,  so  that,  in  principle,  we  may  choose  N  elements  of  the  triangular 
matrices  to  have  any  values  we  please.  In  practice,  it  is  found  that  the 
algorithm  for  determining  the  elements  of  the  triangular  matrices  "grows" 
outwards  from  the  principal  diagonals  of  these  matrices,  and  it  is  convenient 
to  assign  values  to  all  the  elements  on  one  of  the  principal  diagonals; 
usually  they  are  all  made  unity.  In  the  present  application  this  is  not  very 
satisfactory,  because  the  near-symmetry  of  the  matrices  A  and  A^  means  that 
the  triangular  matrices  into  which  A^  is  resolved  will  have  much  the  same 
internal  proportionalities  (the  discussion  in  Reference  2,  pp  142-147,  shows 
that  for  any  symmetric  matrix  there  exists  a  decomposition  into  essentially 
equal  triangular  matrices,  transposes  of  each  other)  but  will  be  markedly 
different  in  their  actual  values,  which  may  lead  to  rounding  errors.  This 
difficulty  has  been  removed  by  adopting  a  triple  decomposition,  setting  every 
element  of  the  principal  diagonals  of  the  triangular  matrices  to  unity  and 
inserting  a  diagonal  matrix  between  the  2  triangular  matrices;  that  is, 

/  lower  triangular  with  \  /  upper  triangular  with 

1  (unit  principal  diagonal/  V  iagona  i  l  n£t  pr£ncipai  diagonal / 
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In  this  decomposition,  the  2  triangular  matrices  are  now  comparable  in  size 
whenever  is  near-symmetric,  they  are  still  banded  and  of  the  same  band¬ 
width  as  Aj  whenever  that  matrix  is  banded,  and  every  element  in  the 
decomposition  that  has  not  already  been  specified  as  unity  or  null  is 
uniquely  determinable  from  the  law  of  matrix  multiplication.  The  uniqueness 
of  this  triple  decomposition  is  proved  in  Appendix  A. 

5  THE  DETERMINATION  OF  THE  ELEMENTS 

We  will  write  the  decomposition  we  intend  to  make,  following  equation  (7), 
in  the  form 

B  =  LDU  (8) 

and  particular  elements  of  these  matrices  will  be  denoted  by  the  correspond¬ 
ing  small  letter  with  suffixes  indicating  the  row  and  column,  eg  b^j  is  the 
element  in  the  i^  row  and  j1"^  column  of  B.  (A^  has  been  replaced  by  B 
merely  to  make  this  convenient  notation  possible.)  We  now  apply  the  law  of 
matrix  multiplication  to  derive  an  expression  for  each  element  in  the  product 
LDU,  which  equation  (8)  demands  shall  be  equal  to  the  corresponding  element 
in  B.  We  have 


(Id)  .  .  =  V  1.,  d,  .  =  1.  .d.  ., 

ij  ik  kj  ij  jj’ 


k  =  l 

since  the  matrices  are  of  order  N  x  N  and  D  is  a  diagonal  matrix:  similarly 


(du)  =  >  d.,u,.=d..u.. 

ij  Z-w  lk  kj  n  ij 


k=l 


Hence  we  may  write 

N  N 

(ldu)  =  N '  (Id)  .u.  =  y '  1  ,d..u. 

rs  A—*  ri  is  t-J  ri  ii  is 

i= 1  i=l 

N  N 

El  . (du) .  -  ^ '  1  .  d .  .  u .  , 

r l  l  s  <  r l  ii  is 

i - 1  i  =  l 

verifying  the  associative  law  of  matrix  multiplication. 
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Consequently  we  have 


k 


t 


* 


1 


t 


1 


N 

b  =  V*  1  .d.  .u. 

rs  n  n  is  (9) 

i~l 

We  now  verify  that  banded  L  and  U  matrices  give  rise  to  a  banded  B  matrix 
of  the  same  bandwidth.  L  is  a  lower  triangular  matrix  with  only  M  non-null 
sub-diagonals,  and  U  is  an  upper  triangular  matrix  with  only  M  non-null 
super-diagonals;  so  we  have 

1  =0  for  r  <  s,  1  =0  for  r  >  s  +  M 

rs  rs 

u  =0  for  r  >  s,  u  =0  for  r  <  s  -  M.  (10) 

rs  rs  '  ' 

From  conditions  (10)  we  obtain  1  .  =  0  for  i  <  r  -  M  and  u.  =0  for  i  >  s; 

n  is 

N 

consequently  all  terms  in  the  sum  Y.  1  .d.  .u.  will  contain  a  null  factor 

.  .  n  n  is 

and  vanish,  making  b  vanish  by  equation  (9) ,  unless  r  -  M  =*  i  .  <  i  =s 

&  rs  j  't  '  min  max 

or  r  <  s  +  M.  Similarly  we  have  1  .  =  0  for  i  >  r  and  u.  =0  for  i  <  s  -  M; 

ri  is 

so  all  terns  in  the  sum  in  equation  (9)  will  contain  a  null  factor  and 

vanish,  making  b  vanish,  unless  s-M=i.  <  i  =  r  or  r  >  s  -  M. 
rs  min  max 

Hence  B  must  have  null  elements  everywhere  except  on  its  principal  diagonal 
and  on  the  M  sub-diagonals  and  M  super-diagonals  closest  to  its  principal 
diagonal;  that  is,  it  must  be  a  banded  matrix  of  bandwidth  M  (or  less). 

This  does  not,  of  course,  prove  that  every  such  banded  matrix  possesses  such 
a  banded  decomposition;  indeed,  no  such  proof  can  exist,  since  it  is  perfectly 
possible  for  a  banded  matrix  to  be  singular  or  to  have  singular  principal 
sub-matrices.  It  can  however  be  proved  that,  if  the  banded  matrix  and  its 
principal  sub-matrices  are  all  non-singular,  a  banded  decomposition  exists, 
which  must  then  be  the  unique  decomposition;  for  the  algorithm  we  are  about 
to  derive  will  always  produce  a  banded  decomposition  provided  that  note  of 
the  diagonal  elements  of  D  vanishes,  and  it  is  shown  in  Appendix  A  (using  a 
proof  based  on  that  of  the  triangular  decomposition  theorem  in  ref  2, 
pp  17-20)  that  this  condition  is  satisfied  whenever  the  banded  matrix  and 
its  principal  sub-matrices  are  all  non-singular. 

It  remains  only  to  obtain  the  unspecified  elements  of  L,  D  and  U,  which 
is  done  by  reversion  of  the  set  of  equations  represented  by  equation  (9). 
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6  THE  REVERSION  ALGORITHM 

We  begin  the  determination  of  an  algorithm  by  expressing  somewhat  more 
compactly  the  information  contained  in  equations  (9)  and  (10)  together  with 
the  conditions  for  the  bandedness  of  B;  thus: 


Known 


Unknown 


i  (r  -  s 
rs 

> 

M)  = 

0 

1  (M 
rs 

>  r 

-  s  >  0) 

-  (r  -  s 
rs 

< 

-  M) 

= 

0 

u  (~ 
rs 

M  < 

r  -  s  <  0) 

(r  -  s 
rs 

> 

M)  = 

0 

d  (r 
rs 

-  s 

=  0) 

(r  -  s 
rs 

< 

0)  = 

0 

(r  -  s 
rs 

= 

0)  = 

1 

(r  -  s 
rs 

> 

0)  = 

0 

rs(r  -  s 

V 

-  M) 

= 

0 

rs(r  "  S 

= 

0)  = 

1 

rs^  -  s 

0)  = 

0 

Determining  Equations 
N 

b  (M  >  r  -  s  >  -  M)  »  V''  1  . d.  .u. 
rs  l-j  ri  n  is 

i=l 

(note  that  the  number  of  determining 
equations  is  the  same  as  the  number  of 
unknown  elements) 

We  introduce  the  two  dyadic  arithmetic  operators  t  and  f ,  which  are  defined 
as  follows: 

x  t  y  «  x  if  x  »•'  y  or  y  if  x  <  y,  ie  the  greater  of  x  and  y 

x  t  y  =  y  if  x  :  y  or  x  if  x  <  y,  ie  the  lesser  of  x  and  y. 

The  summary  just  given  then  becomes,  on  incorporating  the  known  nulls  into 
the  equations  for  the  unknown  elements, 

r-f  s+N 

b  (M  >  r  -  s  >  -  M)  = 

i=  1  t (r-M) t (s-M) 


(1  . d .  . u.  ) 
n  li  is 


(ID 
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(Note  that  the  operators  +  and  +  are  associative  and  comnutative. )  These 

equations  are  most  simply  reversed  by  introducing  the  auxiliary  quantities 

f  defined  by 
rs 

f  =  1  d  ; 
rs  rs  ss 


then  we  obtain,  after  some  simple  algebraic  manipulations. 


f  =  b 
rs  rs 


y]  (f r£u£g)  (s  =l+(r  -  M),  . r  -  ]) 

i=  1  + (r-M) 


d  =  b 
rr  rr 


Y  (f  -u.  ) 
/  n  lr 

i=  1  t (r-M) 


u  =  I  b 
rs  I  rs 


i=  1  t (s-M) 


Zj  (friuis)  /  drr  (s  =  r  +  1 . N+(r  +  M) ) 


1  =  f  /d 

rs  rs  ss 


(+D 

(s  =  1  t(r  -  M) ,  .....  r  -  1) 


b 

where  the  value  of  E  (  ~  )  =  0  whenever  b  <  a.  It  will  be  seen  that 

i=a 

this  algorithm  is  complete  (every  one  of  equations  (11)  is  used,  and  every 

unknown  element  is  calculated),  consistent  (no  element  or  auxiliary  quantity 

is  calculated  more  than  once,  or  referred  to  before  it  has  been  calculated), 

and  will  always  succeed  (the  divisions  in  it  are  always  permissible  because 

none  of  the  diagonal  elements  of  D  can  vanish  under  the  usual  conditions,  and 

all  the  other  operations  are  unconditionally  permissible).  It  can  also  be 

seen  that  the  number  of  arithmetic  operations  required  to  perform  the 
....  2 

decomposition  is  approximately  MN  for  M  <<  N;  this  result,  quoted  in 
Reference  3,  was  confirmed  after  the  algorithm  had  been  programmed,  by 
putting  counting  variables  in  the  loops  of  the  program. 


7  THE  IMPLEMENTATION 

Equations  (8),  (7),  (3)  and  (12)  now  enable  us  to  write 

L]D1U]x(l)  =  f  -  Asx(l_1)  (13) 
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where  the  elements  of  L.. ,  D1 ,  U1  and  A  are  all  known.  If  we  now  determine 

II  I  S  /  •  \ 

by  forward  substitution  the  column  vector  zK  such  that 


L  z(i)  =  f  "  A  x(l-1) 
I  S 


(14) 


(i) 


and  then  determine  by  division  the  column  vector  y  such  that 


_  (i)  _  .  (O 

o xy  =  2 


(15) 


we  then  obtain  by  taking  equations  (13),  (14),  and  (15)  together  the  result 


U^x 


(i) 


=  Y 


(i) 


(16) 


which  may  be  solved  for  by  back-substitution.  The  existence  and 

uniqueness  of  solutions  to  equations  (14),  (15)  and  (16)  is  guaranteed  by 
the  fact  that  the  3  coefficient  matrices  L^,  and  are  non-singular 
whenever  the  sufficient  conditions  for  the  triple  decomposition  of  are 
satisfied. 


From  the  standpoint  of  pure  mathematics,  this  completes  the  problem; 
we  have  specified  an  algorithm  and  determined  the  conditions  under  which  it 
can  be  applied  and  will  converge.  To  complete  the  practical  solution, 
however,  it  is  necessary  to  provide  a  convergence  criterion,  which  is 
usually  based  on  the  value  of  the  residuals  in  the  equations  as  expressed 
by  the  vector  Ax'1  -  f;  from  equations  (1)  and  (4)  this  vector  may  be 
written  as  Ae ^  \  and  we  can  then  use  the  ratio  of  the  magnitudes  of  the 
vectors  Ae^^  and  Ax  (=  f)  as  an  estimate  of  the  ratio  of  the  magnitudes 
ot  the  vectors  e^1^  (=  x^^  -  x)  and  x.  Unfortunately  this  estimate  is  not 
necessarily  a  good  one  if  the  matrix  A  is  ill-conditioned  (ref  2,  pp  103 

et  seq).  Accordingly,  our  implementation  of  the  method  works  by  iterating 

-2 

until  the  error  estimate  is  less  than  some  chosen  small  number  (eg  10  or 
-3 

10  )  and  then  taking  2  more  iterations  and  using  the  3  sets  of  results  so 

obtained  to  obtain  a  final  refined  result,  as  follows.  If  the  error  in  the 
3  sets  of  results  is  determined  solely  by  the  latent  root  largest  in  modulus, 
the  differences  between  them  can  be  manipulated  to  eliminate  the  error 
altogether,  since  the  errors  will  be  decreasing  geometrically.  This 
assumption  is,  of  course,  never  rigorously  true,  and  will  fail  completely 
if  there  happen  to  exist  several  latent  roots  wi th  the  same  largest  modulus 
but  with  different  phases;  but  in  practice  it  is  found  to  give  a  substantial 
increase  in  accuracy,  up  to  an  order  of  magnitude  in  some  cases.  The 


10 


calculation  which  obtains  the  final  result  by  "eliminating"  the  error  also 
provides  us  with  the  value  of  the  convergence  rate,  and  this  value  will  in 
general  be  different  for  each  component  of  the  solution  vector;  hence  the 
consistency  (or  otherwise)  of  the  complete  set  of  convergence  rates  affords 
a  useful  test  of  the  assumption  of  geometric  convergence,  and  if  a 
particular  run  fails  to  converge  we  have  information  from  which  to  estimate 
the  increase  in  bandwidth  required  for  convergence. 

8  PROGRAMMING  DETAILS 

The  algorithm  has  been  implemented  in  Fortran  IV,  following  preliminary 
experiments  mg  Basic,  and  the  Fortran  program  listing  with  explanation  is 
presented  in  Appendix  B.  It  is  assumed  that  the  elements  of  Ag  will  be  held 
in  backing  store,  but  the  elements  of  A^  must  be  kept  in  main  store  because 
the  decomposition  algorithm  involves  reading  A^  both  by  rows  and  by  columns. 
This  restriction  has  prevented  full  testing  of  the  program  to  date,  as 
practical  problems  with  vehicle  antennas  involve  A^  matrices  which  are  too 
large  for  32000  words  of  main  store  (the  maximum  allocated  to  user  programs 
on  many  mainframe  computers,  eg  the  ICL  1900  series  under  George  III); 
however,  tests  with  artificially  constructed  matrices,  and  with  simple  but 
genuine  problems,  have  shown  that  the  method  does  work  and  has  been  implemented 
correctly. 
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APPENDIX  A 


PROOF  OF  THE  TRIPLE  DECOMPOSITION  THEOREM 

It  is  desired  to  prove  that,  if  a  square  banded  matrix  and  all  its 
principal  sub-matrices  are  non-singular,  that  matrix  may  be  written  as  the 
product,  in  order,  of  a  lower-triangular  matrix  of  the  same  bandwidth  and 
with  unit  diagonal  elements,  a  diagonal  matrix  with  no  null  diagonal 
elements,  and  an  upper-triangular  matrix  of  the  same  bandwidth  and  with  unit 
diagonal  elements;  and  that  this  decomposition  is  unique. 

From  reference  2,  pp  17-18,  we  can  see  that  an  arbitrary  square  matrix 
which  is  non-singular  and  has  all  its  principal  sub-matrices  non-singular 
can  be  written  as  the  product,  in  order,  of  a  lower-triangular  matrix  with 
unit  diagonal  elements,  a  diagonal  matrix,  and  an  upper-triangular  matrix 
with  unit  diagonal  elements;  that  this  representation  is  unique;  and  that 
the  diagonal  elements  of  the  diagonal  matrix  are  given  by  the  equations 

lAi  I  lA ! 

d  =  I  A  I  ,  d.  .  (1  <  i  <  N)  -  -  ,  dL^  -  -  , 

11  1  11  lAi-il  iVii 

where  A  is  the  matrix  being  decomposed,  H  is  its  order  (strictly  N  x  N) , 

A.  is  its  itl1  principal  sub-matrix  and  d.  .  is  the  element  in  the  i*"  row 

i  th  11 

and  i  column  of  the  diagonal  matrix  in  the  triple  decomposition  of  A. 

Since  A  and  all  the  A^  are  non-singular,  it  follows  from  these  equations 
that  all  the  d^  are  non-zero.  It  therefore  remains  only  to  prove  that  the 
unique  decomposition  of  a  banded  matrix  A  involves  triangular  matrices  of 
the  same  bandwidth,  and  an  implicit  proof  of  this  is  furnished  by  the 
algorithm  in  the  main  text.  Since  this  provides  a  triangular  decomposition 
(subject  to  the  condition  that  all  the  d^  are  non-zero,  which  we  have  just 
verified),  it  provides  the  (unique)  decomposition,  and  the  decomposition 
involves  banded  triangular  matrices  of  the  appropriate  bandwidth,  which 
completes  the  proof. 
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APPENDIX  B 


THE  FORTRAN  PROGRAM  LISTING,  AND  ITS  EXPLANATION 

The  program  listing  presented  here  divides  into  A  sections:  A  is 
concerned  with  declarations,  initialisation  and  input  of  the  coefficients 
and  righthand  sides;  B  consists  of  the  algorithm  for  decomposing  the  banded 
matrix  into  its  lower-triangular ,  diagonal  and  upper-triangular  components; 

C  contains  the  start  of  the  iteration,  and  its  main  loop;  and  D  contains  the 
iteration  control  code,  the  final  refinement  of  the  results  and  their 
output.  These  sections  will  now  be  explained  in  detail  (we  shall  use  the 
convention  that,  eg  line  A-31  refers  to  the  31st  line  in  section  A). 

In  line  A-l,  the  workspace  arrays  are  declared.  The  sizes  given  allow 
for  a  maximum  matrix  order  of  69  and  a  maximum  bandwidth  of  20  (correspond¬ 
ing  to  A1  non-null  diagonals).  A  is  a  scratch  array,  which  normally  holds 
the  row  of  the  main  matrix  currently  being  operated  on;  B  is  the  array  of 
righthand  sides;  Q  holds  the  non-null  coefficients  of  the  banded  matrix, 
suitably  sheared  in  column  index  so  as  to  convert  a  parallelogram  into  a 
rectangle.  The  variables  QE,  QE1  and  QE2  in  line  A-2  are  scratch  variables, 
as  are  QR,  QI  and  QQ  in  line  A-3;  QR  and  QI  are  generally  used  to  hold  the 
real  and  imaginary  parts  of  QE.  The  variable  SUMBB  in  line  A-3  is  used  to 
hold  the  square  of  the  magnitude  of  the  vector  of  righthand  sides  B,  which 
is  computed  later.  The  arrays  FMTIN  and  FMTOUT  in  line  A-3  are  used  to  hold 
the  object-time  formats  for  input  and  output,  respectively;  the  variable 
ERROR  in  this  line  holds  the  small  number  against  which  the  error  estimate 
is  compared  after  each  iteration.  The  variables  N,  MWIDTH,  ITER  and  MAX 
in  line  A-A  hold  the  order  of  the  matrix  (strictly  this  is  N  x  N) ,  the 
bandwidth,  the  current  iteration  count  and  the  maximum  number  of  iterations 
allowed.  The  variables  I,  J,  K  and  L  in  this  line  are  used  as  indices  for 
loops,  and  the  1A  variables  declared  in  line  A-5  are  scratch  variables  (used 
mostly  for  holding  variable  upper  and  lower  limits  for  loops).  The  variable 
DONE  in  line  A-6  becomes  true  when  the  iteration  is  complete,  and  the 
variable  LARGER  remains  true  as  long  as  the  current  error  estimate  exceeds 
the  value  held  in  ERROR  (as  will  be  explained  in  context,  these  2  conditions 
are  not  complementary). 

In  lines  A- 7  to  A-12,  the  object-time  formats,  the  matrix  order,  the 
permitted  error,  the  permitted  number  of  iterations  and  the  bandwidth  are  read 

B.l 


in,  and  the  tape  or  disc  used  for  out-of-core  storage  is  initialised.  In 
lines  A-13  to  A-26,  the  coefficients  and  the  righthand  sides  are  read  in, 
and  the  coefficients  which  comprise  the  non-null  elements  of  the  banded 
matrix  are  written  to  the  array  Q  (the  indices  are  adjusted  so  that  the 
first  non-null  element  in  each  row  of  the  banded  matrix  is  placed  in  the 
first  element  of  the  corresponding  row  of  Q)  ;  the  coefficients  are  then 
written  out  to  a  scratch  storage  device.  In  lines  A-27  to  A-37,  the  scratch 
storage  file  is  closed  and  reset  (for  the  first  iteration),  the  square  of 
the  magnitude  of  the  vector  of  righthand  sides  is  computed  and  stored  in 
SUMBB,  and  the  vector  of  righthand  sides  B  is  copied  into  the  arrays  X  and 
XB  for  later  use. 

Section  B  consists  essentially  of  one  set  of  nested  loops,  in  which  the 
algorithm  described  in  the  section  "The  reversion  algorithm"  is  executed. 

The  number  of  elements  in  the  matrices  of  the  triple  decomposition  is  the 
same  as  the  number  of  elements  in  the  original  matrix,  if  predetermined  unit 
and  null  elements  are  left  out  of  account,  and  the  algorithm  is  so  arranged 
that  each  element  generated  in  the  decomposition  can  over-write  an  element 
of  the  original  band  matrix  stored  in  Q;  in  most  cases  this  over-writing 
is  performed  more  than  once,  the  intermediate  values  being  used  in  inter¬ 
mediate  calculations.  Each  of  the  non-predetermined  elements  in  the 
triangular  and  diagonal  matrices  is  sheared  in  column  index  so  as  to  preserve 
rectangular  nature  of  the  array  Q;  the  result  of  this  is  that  the  diagonal 
elements  of  the  diagonal  matrix  lie  in  the  (MWIDTH  +  l)t^  column  of  Q,  and 
the  non-predetermined  elements  of  the  triangular  matrices  are  positioned 
touching  on  either  side  (lower-triangular  on  the  left,  upper-triangular 
on  the  right).  With  these  facts  borne  in  mind,  it  can  be  seen  that  section  B 
is  simply  a  translation  of  the  reversion  algorithm  into  Fortran. 

Section  C  begins  by  initialising  several  variables  (lines  C-l  to  05) 

and  then  jumps  from  line  06  to  031,  because  the  first  part  of  the 

iterative  process  is  omitted  on  the  first  iteration  (this  is  equivalent  to 

taking  x^  and  as  null  vectors).  In  lines  07  to  027  we  compute  the 

elements  of  (f  -  and  store  them  in  the  array  X  (lines  09  to  021), 

and  using  the  array  XB  which  contains  (f  -  Agx^  or  (A^x^'  ^)  (see 

equation  (3))  we  compute,  in  lines  022  to  026,  the  square  of  the  magnitude 

of  (Ax^  ^  -  f)  (using  the  consequence  of  equation  (2)  that  (Ax^  ^  -  f)  ■ 

(A.x^-1^)  -  (f  -  A  x"-”))  and  copy  array  X  into  array  XB,  which  is  the 
1  s 
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( i  2) 

reason  why  XB  contains  (f  -  Agx'-  ')  when  it  is  used  the  next  time  round. 

In  lines  C-28  to  C-30  we  compute  the  ratio  of  the  magnitudes  of  the  vectors 
(Ax^  ^  -  f)  and  f,  compare  the  ratio  with  the  permitted  error  held  in 
ERROR,  set  the  logical  variable  LARGER  according  to  the  result  of  the 
comparison,  and  reset  the  file  of  coefficients  for  the  next  iteration. 

The  first  and  subsequent  iterations  now  join  together,  with  array  X 
holding  the  vector  (f  -  A  x^*  ^);  and  the  rest  of  section  C  is  simply  the 
solution  of  equation  (13),  using  the  derived  equations  (14),  (15)  and  (16). 

In  lines  C-31  to  C-41  we  solve  equation  (14)  by  forward  substitution  to  find 
the  vector  ;  in  lines  C-42  to  C-44  we  solve  equation  (15)  by  division 

to  find  the  vector  y^;  and  in  lines  C-45  to  C-56  we  solve  equation  (16) 
by  back-substitution  to  find  the  result  of  the  current  iteration,  the 
vector  x(i).  We  then  update  the  iteration  count  and  output  its  value  and  the 
values  of  2  elements  of  x^,  as  a  check  on  the  progress  of  the  iteration. 

Lines  D-l  and  D-2  contain  tests  for  the  completion  of  the  iteration. 

Normally  both  these  tests  will  fail,  and  control  will  pass  to  line  D-3;  in 

lines  D-3  to  D-6,  the  solution  just  computed  will  be  copied  into  the  array 

reserved  for  the  previous  solution  and  control  will  be  returned  to  line  C-7 

where  the  next  iteration  begins.  However,  if  on  the  previous  iteration  (not 

the  one  just  completed  but  the  one  before  that)  the  error  estimate  has 

achieved  or  improved  on  the  permitted  error,  the  test  in  line  D-2  will 

succeed  and  control  will  pass  to  line  D-7;  this  will  also  happen  if  the 

iteration  just  completed  has  brought  the  iteration  count  up  to  the  maximum 

allowed  number.  In  lines  D-7  to  D-12,  the  previous  solution  and  the  new 

solution  are  both  moved  one  level  to  become  the  previous-but-one  solution 

and  the  previous  solution  (respectively),  the  logic  variable  DONE  is  set  to 

indicate  that  the  coming  iteration  is  to  be  the  last  one,  and  control  is 

returned  to  line  C-7  for  this  last  iteration.  When  it  is  complete,  control 

arrives  at  line  D-l  as  usual,  and  this  time  the  test  in  this  line  is 

successful  and  control  passes  to  line  D-13.  Assuming  that  by  this  stage  the 

error  in  each  of  the  unknowns  is  determined  wholly  or  mainly  by  the  latent 

root  largest  in  modulus  of  the  iteration  matrix  (-  A.  ^A  ) ,  we  may  write 

X  s 

corresponding  elements  of  the  3  solutions  we  now  have  as  (a  +  6)  (previous- 

2 

but-ono),  (a  +  gX)  previous)  and  (a  +  BA  )  (present),  where  a  is  the  true 

value  and  A  is  the  latent  root  largest  in  modulus.  In  lines  D-15  and  D-16 
2 

the  values  of  (gX  -  gX)  and  (gX  -  g)  are  computed  and  assigned  to  the 


variables  QE1  and  QE2  respectively.  If  (BA  -  B)  is  found  to  be  zero  within 

the  resolution  of  the  computer,  indicating  that  6  is  essentially  zero  or  A 

is  essentially  unity,  the  appropriate  element  of  array  A  is  zeroed  (as  a 

marker)  and  the  present  solution  in  array  X  is  left  unchanged  (lines  D-17  to 

D-19);  this  is  consistent  with  the  fact  that  if  tf  «  0  or  A  *  1  the  sequence 

2 

of  solutions  (i  +  l),  (a  +  BA),  («  +  BA  )  cannot  be  improved  on.  If 
(f\  -  8)  is  not  essentially  zero,  control  passes  from  line  D-17  to  D-20,  and 
in  lines  D-20  and  D-21  we  then  compute  A  and  u  and  assign  them  to  the 
appropriate  elements  of  arrays  A  and  X,  using  the  identities  A  ■■ 

(8>2  -  BA)/(PX  -  B)  and  a  -  (a  +  FA2)  -  ((BA2  -  BA)2/((BA2  -  BA)  -  (BA  -  t))). 
Finally,  in  lines  D-23  to  D-26  we  write  out  the  results  and  terminate  the 
program  run. 


t 


* 


a 
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SUPPLEMENT  TO  APPENDIX  B  THE  FORTRAN  PROGRAM  LISTING 


PART  A 


1 .  COMPLEX  A  (69)  »  E<(69)  »X(6?)  » XL  AST  (69)  »XLAST2(69)  » XB( 69 ) » G ( 69 »  41 ) 

2.  COMPLEX  GE.QE1 .QE2 

3.  REAL  FMTIN(8) ,FMT0UT(8> , ERROR » SUMBB » QR , QI , QQ 

4.  INTEGER  N » MWIDTH » ITER » MAX » I » J » K  »L 

5.  INTEGER  1 1 »  J1 »  J2 . J3» K1 » K2 » Ml » M2. M3. IM1 » JM1 . KM.KM1 » KM2 

6.  LOGICAL  DONE. LARGER 

7.  READ  (500.1000)  FMTIN 

8.  READ  (500.1000)  FMTQUT 

?.  1000  FORMAT  ((6X.8A8)) 

10.  READ  (500.1010)  N, ERROR. MAX. MWIDTH 

11.  1010  FORMAT  ( 14. F10. 0.214) 

12.  REWIND  25 

13.  DO  20  1=1. N 

14.  READ  (550)  ( A  <  L  >  »  L  =  1.N> 

15.  READ  (500. FMTIN)  ( B  ( L. ) .  L  =  1.N) 

16.  J1=I -MWIDTH 

17.  t  M 1  =  1  J 1 

18.  IF  (JJ  .LT.  1)  Jl  =  l 

19.  J2= I TMWI DTH 

20.  IF  ( J2  .GT.  N)  J2=N 

21.  DO  10  J=J1.J2 

22.  JM1 = IM1 f J 

23.  Q(  I ,  JM1  )=A(  .1) 

24.  10  CONTINUE 

25.  WRITE  (25)  ( A ( L  > .  L  =  1.N) 

26.  20  CONTINUE 

27.  ENDFILE  25 

28.  REWIND  25 

29.  SUMBB =0.0 

30.  DO  30  1=1. N 

31.  GE=B  ( I ) 

32.  GR=REAL (GE ) 

33.  QI=AIMAG( GE ) 

34.  SUMBB=SUMBB+GR*GRfGI«GI  . 

35.  X ( I ) =QE 

36.  XB ( I )  =QE 

37.  30  CONTINUE 


SUPPLEMENT  TO  APPENDIX  B  -  THE  FORTRAN  PROGRAM  LISTING 


PART  B 


1 . 

M1=MWIDTH+1 

2. 

DO  170  1=1 

3. 

J1--I  -MWIDTH 

4. 

IF  <J1  .LT.  1)  Jl=i 

5  * 

IM1=M1-  I 

6  ♦ 

DO  110  J=J1.I 

7. 

K2=J-1 

8. 

IF  <J1  .G T.  K2)  GOIO  110 

9. 

JM1=IM1+J 

10. 

KM=M1+J 

11 . 

QE=Q( I.JM1) 

12. 

DO  100  K=J1.N2 

13. 

KMl  =  IMHK 

14. 

KM2=KM  N 

15. 

QE-'GE-G  ( I » KM1 )»Q<N.NM2> 

16. 

100 

CONTINUE 

17. 

0(1 . JM1 )=QE 

18. 

110 

CONTINUE 

19. 

J3=I 1MWIDTH 

20. 

IF  (J3  . GT .  N)  J3=N 

21 . 

J2=I f 1 

■i? 

K2=  I "  1 

23. 

IF  <J2  .GT.  J3>  GOTO  150 

24. 

DO  140  J~J2> J3 

25. 

K1=J- MWIDTH 

26. 

IF  (K1  .LT.  1)  Nl  =  3 

27. 

JM1=IM1+J 

28. 

KM=M1+J 

29. 

QE  =  0 ( I  *  JM1 ) 

30. 

IF  ( N 1  .GT.  K2  >  GOTO  130 

31 . 

DO  120  K=N1 *  N2 

32. 

NM1  =IMH  N 

33. 

NM2=KM- N 

34. 

GE ■ QE- G  < I f  NM1 ) #Q  <  K » NM7 ) 

35. 

120 

CONTINUE 

36. 

130 

Q( If JM1 ) =  QE/G ( I » Ml ) 

37. 

140 

CONTINUE 

38. 

150 

IF  iJl  .GT.  N2 )  GOTO  170 

39. 

DO  160  J  =  J 1 » K 2 

40  . 

IM1  =  IM1+J 

41  . 

CM  I f JM1 )=Q( I • JM1 >/G< J»M1 > 

42. 

160 

CONTINUE 

43. 

170 

CONI INUE 

r 


J 
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1 . 

I TER=0 

n 

aL  4 

LARGER=, TRUE. 

3. 

DONE =. FALSE . 

4  . 

M2=M1+1 

c r 

J  ♦ 

M3=M1 FMWIDTH 

6  * 

GOTO  260 

-t 

t  * 

200 

00=0.0 

0. 

DO  250  1=1, N 

9. 

Q£=B ( I ) 

10, 

READ  (25)  <  A ( L ) »  L=1»N> 

11 . 

J1=I~M1 

12. 

IF  <J1  .LT.  1)  GOTO  220 

13. 

DO  210  J=1 , J1 

14. 

QE  =  QE  A ( J ) #XLAST  <  J ) 

15. 

210 

CONTINUE 

16. 

220 

J2=I EMI 

17. 

IF  ( J2  .GT.  N)  GOTO  240 

18. 

DO  230  J=J2,N 

1?. 

QE=QE  A<  J)*XL.AST<  J) 

20. 

230 

CONTINUE 

21  . 

240 

X  < I ) =  QE 

a:.  • 

OE=XB ( I )  OE 

23. 

QR=REAL ( Of  > 

24. 

QI =AI MAG ( QE ) 

25. 

OQ=OQEQR*OR+OT  #0  I 

26. 

XB< I >=X( I > 

27. 

250 

CONTINUE 

28. 

00= SORT ( QO/SUM/D) 

29. 

LARGER  *  00  .GT,  ERROR- 

30. 

REWIND  25 

31 , 

260 

DO  280  1  =  2  >i, 

32. 

OE  =  X< I  ) 

33. 

J1 =M2~ I 

34  . 

ir  <ji  .Li.  n  ji=i 

35. 

IM1=I  Ml 

36. 

DO  2  70  I--  Jl  rMWIDTH 

37. 

JM 1 = I M 1 f J 

38. 

QE  =  QE  0(I,.D*X(  JM1  ) 

3? . 

2  70 

CONTINUE 

40. 

X ( I >=0E 

41  . 

200 

CONTINUE 

42. 

DO  290  I  1  ,N 

43. 

<  <  I  >  X  (  1  >  /  0  t  J  •  M  1  ) 

44  . 

290 

CONTINUE 

45. 

DO  310  ]  “  .  ‘ ,  N 

46. 

1 1  =n  1 1  r 

47. 

GC  =  X< 11  ) 

40. 

J3=I+MWID  III 

49. 

IF  <J3  . G 7 .  M3  >  J3=M3 

50. 

I Ml = I 1 -Ml 

51  . 

DO  300  J=M2 » J3 

r  "> 

J*.  » 

JM 1  =  I M 1 f  J 

j3. 

QE=Or  0(11 f J>*X< JM1 ) 

54. 

300 

CONTINUE 

55  . 

1(11 >  =GE 

56. 

310 

CONTINUE 

57. 

ITER-  ITER)  J 

58. 

WRITE  ( 600, FMTOUT)  [TER 
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1,  IF  (DONE)  GOTO  430 

2,  IF  ((.NOT.  LARGER)  .OR.  (ITER  .EG.  MAX))  GOTO  410 

3.  DO  400  1  =  1,  N 

4.  XLAST ( I ) =X( I ) 


5 . 

400 

CONTINUE 

6 . 

GOTO  200 

7. 

410 

DO  420  1  =  1,  N 

8. 

XLAST2( I )=XLAST (  I ) 

9. 

XLAST 1 1 >  =X ( I ) 

10. 

420 

CONTINUE 

11 . 

DONE =  . TRUE. 

12. 

GOTO  200 

13. 

430 

DO  460  I  =  1,N 

14. 

QE=X ( I ) 

15. 

GEJ.  =GE- XLAST  ( I ) 

16. 

QE2=GE~QE1  -XL.AST2  ( I ) 

17. 

IF  ( CADS ( GE2 ) >  450,440,450 

18. 

440 

A  ( I )  =  (  0 . 0  ,  0 . 0 ) 

19. 

GOTO  460 

20. 

450 

A ( I ) =QE 1 /QE2 

21 . 

X(I)  =  GE  -  QE1#<QE1/(QE1-GE2> ) 

ti.  4 

460 

CONTINUE 

23. 

WRITE  (600,1020)  MWIDTH , I TER , MAX , ERROR , GG 

24 , 

1020 

FORMAT  ( 1 H 1 >  3 1 5 , 2X  >  2 ( 2X , E 1 4 . 7 ) ) 

25. 

WRITE  ( 600 , FMTOUT )  ( ( L  *  A ( L  > , X ( L ) ) »  L  =  1»N> 

26 . 

STOP 

27. 

END 
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